home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / VideoToolbox 96.06.15 / VideoToolboxSources / Luminance.c < prev    next >
Text File  |  1996-03-07  |  41KB  |  1,077 lines
  1. /*
  2. Luminance.c
  3. See Luminance.h for prototypes.
  4. See Luminance.note for explanation.
  5. Also see:
  6. D.G. Pelli and L. Zhang (1991) Accurate control of contrast on microcomputer displays. 
  7. Vision Research, 31:1337-1360.
  8.  
  9. Copyright (c) 1989-1995 Denis G. Pelli. This software is free; you may use it in
  10. your research and give it away to others, with the following restrictions. Any
  11. copy you give away must include this paragraph, unmodified, and, if you have
  12. made any changes to the rest of the file, must include a note, added to HISTORY
  13. below, giving your name, the date, and a description of the changes. This
  14. software may not be sold, whether in source or compiled form, without my
  15. permission. I hope you will find this software useful, but I can't promise that
  16. it will work for you, and am not offering any support. That's why it's free. I
  17. would appreciate reports of bugs and improvements.
  18.  
  19. Denis G. Pelli
  20. Department of Psychology
  21. New York University
  22. 6 Washington Place
  23. New York, NY 10003
  24. denis@cns.nyu.edu
  25.  
  26. HISTORY:
  27.  
  28. 4/24/89 dgp first working version.
  29. 4/28/89 dgp added device argument.
  30. 6/23/89 dgp added support for a video attenuator that has unequal gain in the three
  31.             channels. Added LoadClut() and LuminanceClut() routines.
  32. 7/30/89 dgp replaced call to SetEntries (which calls the Color Manager) by a call
  33.             to GDSetEntries (which calls the device driver).
  34. 8/10/89 dgp fixed bugs.
  35. 8/13/89 dgp rewrote LToV() to use nth order polynomial. Broke up LinearizeClut into
  36.             three routines--SetLuminanceRange, SetLuminance, and SetLuminances--the
  37.             last of which is equivalent to the old LinearizeClut.
  38. 8/16/89    dgp    After reading the Brooktree DAC manual I changed the documentation and
  39.             increased tolerance from 0.5 to 1 LSB, and now use the highest
  40.             luminance gain (worst case) over the requested range to transform it
  41.             into a luminance difference tolerance.
  42. 8/21/89    dgp    Made minor improvements in LToV() to deal with failure to converge due to 
  43.             getting stuck in a local minimum of the quartic. My solution is to now
  44.             require that Lmin and Lmax be filled in by the user in the calibration
  45.             structure. In practice this will guarantee that the minimum luminance
  46.             will not be below the local minimum found at the base of the rising curve.
  47.             Last week I also changed the LToV algorithm slightly, to do a bisection if
  48.             Newton's method would take us out of range.
  49. 9/7/89 dgp    Fixed bug in sorting routine, Sort().
  50. 9/10/89 dgp    Commented out the polynomial versions of VToL and LToV and replace them by
  51.             new versions that use a power law. The power law is a marginally better fit
  52.             than a 6th order polynomial (using only 4 instead of 7 parameters) and its
  53.             inverse can be computed much more quickly, about 0.3 ms instead of 2 ms.
  54. 9/11/89 dgp    Deciding to be indecisive, I introduced a conditional POWER_LAW_FIT
  55.             to control the choice of power law or polynomial fit for LToV and VToL.
  56. 9/16/89    dgp    Having finished writing the paper (Pelli & Zhang, 1991) I came back to
  57.             change the program to agree with what I wrote. I changed the equivalent
  58.             number tolerance to be the sum instead of the maximum of the 
  59.             variable-dac gains.
  60. 9/25/89 dgp Fixed minor bug in computation of tolerance. Now gives correct answer even
  61.             when the lowLuminance or highLuminance is out of range.
  62. 10/9/89 dgp    Made minor change so that when some of the gain[] are zero SetLuminance will
  63.             load zero into the unused colorSpec array elements. This is only a cosmetic
  64.             change.
  65. 10/28/89 dgp Made minor changes. I declared LToVPolynomial and LToVPower and
  66.             VToLPolynomial and VToLPower instead of having a conditional compilation. 
  67.             This makes both flavors of routine always available. I also changed 
  68.             LToVPolynomial to use LToVPower instead of LToVQuadratic for its initial 
  69.             guess. The quadratic guess was often poor and occasionally led to 
  70.             convergence problems. If the power law fit is not available then it reverts
  71.             to using the quadratic fit. If that's not available then it reverts to using
  72.             a middle-of-range guess.
  73. 10/28/89 dgp I eliminated the compile-time flag POWER_LAW_FIT and instead use whichever, 
  74.             polynomial or power, provides a better fit, as determined at run time. I 
  75.             biased the comparison of rms error to favor the power law fit since it has 
  76.             few parameters and can be inverted much more quickly. Note: if the 
  77.             power, polynomial, or quadratic fit parameters are not supplied then the 
  78.             appropriate LR.powerError, LR.polynomialError, or LR.quadraticError field 
  79.             should be set to infinity: INF or 1.0/0.0.
  80. 12/5/89    dgp    Added the trivial routine LToL() which enforces the bounds LP->LMin and 
  81.             LP->LMax.
  82. 4/2/90    dgp    Capitalized the word "To" in all function names, e.g. LToL().
  83. 4/22/90    dgp    Version 1.4. Made minor changes to the documentation. Renamed LToV() to
  84.             LToVFormulaic(), and introduced a new LToV() that uses a table to run faster.
  85.             LToV() takes an average of 170 µs, whereas LToVFormulaic() takes
  86.             1700 µs. I also now use the LuminanceTable, if available, to speed up
  87.             VToL() slightly, from 310 µs to 180 µs. Each call to SetLuminance() now
  88.             takes 0.9 ms whereas it used to take 2.5 ms. The loss in accuracy is
  89.             negligible. The interpolation error can be reduced still further by
  90.             increasing LUMINANCES_IN_TABLE. Each doubling of the table size quarters the
  91.             maximum possible error of the interpolation.
  92. 4/23/90    dgp    Added a few lines of code to LToV() to check near the last index, in the
  93.             spirit of Numerical Recipes in C, hunt.c, before starting the bisection
  94.             search. I updated the times in the 4/22/90 note to reflect the latest timing
  95.             by TestLuminance.
  96. 7/27/90    dgp    Tightened up the error checking for illegal entries into the clut. I have
  97.             the impression that, contrary to specification, the Apple video driver
  98.             crashes and corrupts itself if given an out-of-range entry value in a
  99.             setEntry Control call.
  100. 9/18/90    dgp    Changed all instances of "v" to "V". The final version of the Pelli
  101.             & Zhang (1991) paper refers to the nominal voltage v; this file
  102.             now refers to the "equivalent number" V; they are related by V=255*v.
  103. 9/24/90    dgp    Updated the documentation to correspond more closely to the (hopefully)
  104.             final version of the Pelli & Zhang (1991) paper.
  105. 10/29/90 dgp Doubled speed of SetLuminance(). Now SetLuminance() takes 133 ms to
  106.             do a complete clut for a small new luminance range, and takes 188 ms for a
  107.             large new range. This required minor changes to SetLuminanceRange() & LToV().
  108.             Changed all function headers to Standard C prototype style.
  109. 10/30/90 dgp Replaced phrase "clut value" by "nominal voltage", for consistency with
  110.             Pelli & Zhang (1991).
  111.             Introduced new fixed fraction data type called "Milli", defined in Milli.h.
  112.             By judicious replacement of double by Milli, particularly in the luminance
  113.             table, I have increased the speed of SetLuminance by a factor of three.
  114.             A Mac II now takes only 42 to 50 ms to build a whole clut. This new feature
  115.             is enabled by setting FAST_LUMINANCE to 1 in Luminance.h. There is a slight
  116.             increase in the luminance tolerance.
  117. 10/31/90 dgp More fine tuning. Now takes 31 to 44 ms to build a whole clut. Introduced
  118.             conditional FAST_LUMINANCE in Luminance.h that causes the same code to
  119.             be compiled with either Milli or double variables. The _SetLuminance()
  120.             routine is now quite polished, and runs fast. All variables and routines
  121.             whose names begin with underscore _ are written here with type Milli, but
  122.             if FAST_LUMINANCE is false then this type is redefined (in this file only)
  123.             as double, and all the Milli macros are redefined to be appropriate for
  124.             a double. So bear in mind that things beginning with underscore are of
  125.             unknown type. The virtue of being able to run the same code with Milli or
  126.             double computations is that if you suspect an error in the (homemade)
  127.             Milli arithmetic, you can try out double arithmetic simply by setting
  128.             FAST_LUMINANCE to 0 and recompiling.
  129.             In the interest of speed I have streamlined the tolerance
  130.             computation. The new tolerances are probably as good as the old ones.
  131. 11/2/90 dgp Added _SetLuminances() subroutine that substitutes linear interpolation
  132.             for most of the calls to _LToV(). This has greatly speeded up
  133.             SetLuminancesAndRange(), with minor loss of accuracy. A complete clut now
  134.             takes 9 to 31 ms. Even a Mac II can now compute low-contrast cluts
  135.             on the fly, between frames. (The user may wish to optimize the speed-
  136.             accuracy trade-off controlled by the LINEAR_V_DOMAIN parameter, which
  137.             determines the maximum width of the interpolation interval.)
  138. 11/6/90 dgp Replaced Milli by FIXED, which can be compiled as either double or Fixed.
  139.             The Fixed math is sufficiently precise as to give the same DAC values as
  140.             double math.
  141.             Now figure out optimal bit shift LT->L.LShift to preserve as much resolution
  142.             as possible when interpolating in _LToV(). 
  143.             A complete clut now takes 8 to 30 ms. For threshold stimuli the average
  144.             time will usually be less than 15 ms.
  145. 11/8/90 dgp    Eliminated gamma slope table since the speed-up it offered was too small
  146.             to measure. Tidied up the computation of LShift.
  147. 8/24/91    dgp    Made compatible with THINK C 5.0.
  148. 12/17/91 dgp Added new routines: IncrementLuminance() and GetV().
  149.             The former is used to obtain the lowest possible contrast. On a log contrast
  150.             scale, minus infinity (zero contrast) is a poor approximation to even
  151.             a small log contrast.
  152. 3/11/92    dgp    Minor editing of comments.
  153. 12/2/92 dgp Fixed minor bug in SetLuminanceRange() reported by Wei Xie and Ken Alexander
  154.             that could cause the THINK C Debugger to report an error when an 
  155.             uninitialized double V[] was converted to an int. The value was not used 
  156.             for anything, so this was innocuous.
  157. 12/17/92 dgp 1.94 Began enhancement to allow any dacSize, not just 8-bit dacs. I have
  158.             confirmed that everything still works fine with 8-bit dacs, but I don't
  159.             have a 9-bit dac video card to test it any further. •CAUTION: this probably
  160.             does not yet work correctly for 9-bit dacs. E.g. it seems likely that the 
  161.             fixed point math will break if VMax is increased to 511; a lot of careful 
  162.             error analysis went into the original coding, at which time I thought that 
  163.             it was safe to assume an 8-bit dac (i.e. VMax==255). •Removed obsolete 
  164.             support for THINK C 4. 
  165. 12/21/92 dgp No longer load unused dac bits.
  166. 1/4/93    dgp    LoadLuminances now checks the gdType.
  167. 5/10/93    dgp    If driver is in gray mode (i.e. not color) then LoadLuminances maps rgb to 
  168.             gray by simply copying the green component to the red and blue components.
  169. 5/12/93    dgp removed debugging code that forced isGray to always be true.
  170. 3/24/94    dgp    added LToE, LToEOrdered, and EToL.
  171. 6/29/94 dgp renamed "LtoEOrdered" to "LToEOrdered", fixing erroneous capitalization.
  172. 9/28/94    dgp    If necessary, IncrementLuminance now bumps up the upper luminance bound to 
  173.             include the incremented luminance.
  174. 10/25/94 dgp LoadLuminance and GetV now ignore a device with no driver. (E.g. a device
  175.             created by NewGWorld.)
  176. 12/29/94 dgp Fixed zero-divide halt in computing FAST_LUMINANCE _LToV returned value,
  177.             reported by Josh Solomon. Increased JMAX from 20 to 30 to eliminate spurious
  178.             warning in LToVPolynomial.
  179. 3/5/96 dgp Fixed error in GetV(), which called GDGetEntries() incorrectly.
  180. */
  181. #include "VideoToolbox.h"        /* prototype for GDSetEntries() */
  182. #include "Luminance.h"
  183. static void Sort(int n,double arr[],int krr[]);    /* for internal use only */
  184.  
  185. #undef LongToFix
  186. #undef FixToLong
  187. #undef DoubleToFix
  188. #undef FixToDouble
  189. #if FAST_LUMINANCE
  190.     #define FIXED Fixed
  191.     #if !MACINTOSH
  192.         typedef long Fixed;
  193.         #define FixMul(x,y) DoubleToFix(FixToDouble(x)*FixToDouble(y))
  194.         #define FixDiv(x,y) DoubleToFix(FixToDouble(x)/FixToDouble(y))
  195.     #endif
  196.     #define LongToFix(x) ((long)(x)<<16)
  197.     #define FixToLong(x) ((x)>>16)
  198.     #define DoubleToFix(x) ((Fixed)((x)*65536.+0.5))
  199.     #define FixToDouble(x) ((double)(x)*(1./65536.))
  200.     #define D(x) FixToDouble(x)                            /* handy for debugging */
  201. #else
  202.     #define FIXED double
  203.     #define FixMul(x,y) ((double)(x)*(y))
  204.     #define FixDiv(x,y) ((double)(x)/(y))
  205.     #define LongToFix(x) ((double)(x))
  206.     #define FixToLong(x) ((long)(x))
  207.     #define DoubleToFix(x) ((double)(x))
  208.     #define FixToDouble(x) ((double)(x))
  209. #endif
  210.  
  211. #undef MAX
  212. #undef MIN
  213. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  214. #define MIN(a,b) ((a) < (b) ? (a) : (b))
  215.  
  216. double EToL(LuminanceRecord *LP,int entry)
  217. {
  218.     return GetLuminance(NULL,LP,entry);
  219. }
  220.  
  221. int LToE(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  222. //Returns the index of the table entry in specified range with luminance closest to L.
  223. {
  224.     int i,entry;
  225.     double dL,latestDL;
  226.     
  227.     latestDL=INF;
  228.     entry=0;
  229.     for(i=firstEntry;i<=lastEntry;i++){
  230.         dL=fabs(L-GetLuminance(NULL,LP,i));
  231.         if(dL<latestDL){
  232.             entry=i;
  233.             latestDL=dL;
  234.         }
  235.     }
  236.     return entry;
  237. }
  238.  
  239. int LToEOrdered(LuminanceRecord *LP,double L,int firstEntry,int lastEntry)
  240. /*
  241. Returns the index of the table entry in specified range with luminance closest to L.
  242. Assumes ordered table from firstEntry to lastEntry, either increasing or decreasing.
  243. Returned entry is guaranteed to be in the range firstEntry...lastEntry.
  244. Based on Numerical Recipes "locate.c", page 117
  245. */
  246. {
  247.     register int lower,middle,upper;
  248.     register double LLower,LMiddle,LUpper;
  249.     int ascend;
  250.     
  251.     lower=firstEntry;
  252.     upper=lastEntry;
  253.     LLower=GetLuminance(NULL,LP,lower);
  254.     LUpper=GetLuminance(NULL,LP,upper);
  255.     ascend=(LUpper>LLower);
  256.     while(upper-lower>1){
  257.         middle=(upper+lower)/2;
  258.         LMiddle=GetLuminance(NULL,LP,middle);
  259.         if(L>LMiddle == ascend){
  260.             lower=middle;
  261.             LLower=LMiddle;
  262.         }
  263.         else{
  264.             upper=middle;
  265.             LUpper=LMiddle;
  266.         }
  267.     }
  268.     if(L-LLower<LUpper-L == ascend)return lower;
  269.     else return upper;
  270. }
  271. double SetLuminance(GDHandle device,LuminanceRecord *LP
  272.     ,int theEntry,double luminance
  273.     ,double lowLuminance,double highLuminance)
  274. /*
  275. Set one entry in the ColorSpec table (and the clut if device is not NULL) to
  276. a specified luminance. It's ok for lowLuminance to be greater than highLuminance; 
  277. they still designate a range.
  278. */
  279. {
  280.     FIXED _luminance;
  281.     
  282.     if(theEntry<0 || theEntry>=COLORS 
  283.         || device!=NULL && theEntry>=GDCLUTSIZE(device))
  284.             PrintfExit("\007SetLuminance: illegal entry %d\n",theEntry);
  285.     
  286.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  287.     _luminance=DoubleToFix(luminance);
  288.     _SetLuminance(LP,theEntry,_luminance);
  289.     if(device != NULL)LoadLuminances(device,LP,theEntry,theEntry);
  290.     return FixToDouble(_Tolerance(LP,_luminance));
  291. }
  292.  
  293. double SetLuminances(GDHandle device,LuminanceRecord *LP
  294.     ,int firstEntry,int lastEntry
  295.     ,double firstLuminance,double lastLuminance)
  296. /*
  297. Set a series of entries in the ColorSpec table (and the clut if device is not NULL)
  298. to a linear sequence of luminances. Assume this is the entire luminance range of 
  299. interest.
  300. */
  301. {
  302.     return SetLuminancesAndRange(device,LP,firstEntry,lastEntry
  303.         ,firstLuminance,lastLuminance,firstLuminance,lastLuminance);
  304. }
  305.  
  306. double SetLuminancesAndRange(GDHandle device,LuminanceRecord *LP
  307.     ,int firstEntry,int lastEntry
  308.     ,double firstLuminance,double lastLuminance
  309.     ,double lowLuminance,double highLuminance)
  310. /*
  311. Set a series of entries in the ColorSpec table (and the clut if device is not
  312. NULL) to a linear sequence of luminances. Uses last two arguments to set the
  313. luminance range.
  314.  
  315. I introduced a stack-space check before calling _SetLuminances(), since it calls
  316. itself recursively and may use 1000 bytes all together. However, printing an
  317. error message requires at least 4500 bytes. So I quit if the stack gets below
  318. 6000, so that the user will be presented with an understandable error message
  319. rather than a mysterious quit or crash.
  320. */
  321. {
  322.     FIXED _tolerance,_t,_dL64;
  323.     FIXED _firstL,_lastL,_firstV,_lastV;
  324.     
  325.     if(firstEntry<0 || firstEntry>lastEntry || lastEntry>=COLORS 
  326.         || device!=NULL && lastEntry>=GDCLUTSIZE(device) )
  327.         PrintfExit("\007SetLuminancesAndRange: illegal entries %d %d\n",firstEntry,lastEntry);
  328.     if(StackSpace()<6000)PrintfExit("SetLuminancesAndRange: not enough stack space.\007\n");
  329.     SetLuminanceRange(LP,lowLuminance,highLuminance);
  330.     _firstL=DoubleToFix(firstLuminance);
  331.     _lastL=DoubleToFix(lastLuminance);
  332.     _firstV=_LToV(LP,_firstL + LP->_LOffset);
  333.     _lastV=_LToV(LP,_lastL + LP->_LOffset);
  334.     if(lastEntry != firstEntry){
  335.         #if FAST_LUMINANCE
  336.         _dL64=DoubleToFix(64.0*(lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  337.         #else
  338.         _dL64=DoubleToFix((lastLuminance-firstLuminance)/(lastEntry-firstEntry));
  339.         #endif
  340.     }else _dL64=0;
  341.     _SetLuminances(LP,firstEntry,lastEntry,_firstL,_dL64,_firstV,_lastV);
  342.     if(device != NULL)LoadLuminances(device,LP,firstEntry,lastEntry);    
  343.     
  344.     /* quick and dirty estimate of tolerance */
  345.     _tolerance=LP->L._Lu[0] - _firstL;
  346.     _t=_lastL - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  347.     if(_tolerance<_t)_tolerance=_t;
  348.     if(_tolerance<0)_tolerance=0;
  349.     _tolerance+=LP->_tolerance;
  350.     return FixToDouble(_tolerance);    /* estimate of max error in ΔL */
  351. }
  352.  
  353. void _SetLuminances(LuminanceRecord *LPtr,int first,int last
  354.     ,FIXED _firstL,FIXED _dL64,FIXED _firstV,FIXED _lastV)
  355. /*
  356. This routine eliminates most of the the slow _LToV() table lookups and
  357. interpolations by assuming that the gamma function is smooth enough that we may
  358. linearly interpolate over any interval of V no wider than LINEAR_V_DOMAIN. Since
  359. the gamma function is roughly parabolic, the consequent error in L will be
  360. proportional to the square of LINEAR_V_DOMAIN, and independent of where this
  361. interval is along V. Since we're inscribing straight line segments in a
  362. positively accelerating gamma function, we will overestimate luminance, and
  363. consequently V will be too low. The error will be greatest at the middle of each
  364. linearly interpolated V interval.
  365.  
  366. If the requested V interval is larger than LINEAR_V_DOMAIN, then it is split
  367. into two intervals by making a pair of recursive calls. LINEAR_V_DOMAIN is
  368. defined in Luminance.h. I suggest a value of 4, but it may be set larger to
  369. attain slightly higher speed at lower accuracy.
  370.  
  371. CAUTION: This routine is a stack hog. Each call uses up about 64 bytes on the
  372. stack, and it calls itself recursively, up to log2(256/LINEAR_V_DOMAIN) times.
  373. Do NOT change any of the declarations of the "new" variables used in the first
  374. if{} to "static", as that would cause the recursion to fail.
  375. */
  376. {
  377.     static int doLast=1;
  378.     int saveDoLast;
  379.     int newOne;
  380.     FIXED _newL,_newV;
  381.     register LuminanceRecord *LP=LPtr;
  382.     register int i,Vi;
  383.     register FIXED _VToGo,_g;
  384.     static RGBColor *RGBPtr;    /* static in order to minimize stack usage */
  385.     static int theEntry;        /* static in order to minimize stack usage */
  386.     static FIXED _V,_dV;        /* static in order to minimize stack usage */
  387.     register short leftShift;
  388.     
  389.     if(last-first>1 ){
  390.         _dV=_lastV-_firstV;
  391.         if(_dV>LongToFix(LINEAR_V_DOMAIN) || _dV<LongToFix(-LINEAR_V_DOMAIN)){
  392.             newOne=(first+last)>>1;
  393.             #if FAST_LUMINANCE
  394.                 _newL=_firstL+((newOne-first)*_dL64>>6);
  395.             #else
  396.                 _newL=_firstL+(newOne-first)*_dL64;
  397.             #endif
  398.             _newV=_LToV(LP,_newL + LP->_LOffset);
  399.             saveDoLast=doLast;
  400.             doLast=0;
  401.             _SetLuminances(LP,first,newOne,_firstL,_dL64,_firstV,_newV);
  402.             doLast=saveDoLast;
  403.             _SetLuminances(LP,newOne,last,_newL,_dL64,_newV,_lastV);
  404.             return;
  405.         }
  406.     }
  407.     if(last!=first)_dV=(_lastV-_firstV)/(last-first);
  408.     else _dV=0;
  409.     if(!doLast)last--;
  410.     leftShift=LP->leftShift;
  411.     for(theEntry=first,_V=_firstV;theEntry<=last;theEntry++,_V+=_dV){
  412.         /****** This section of code is copied from _SetLuminance() ****/
  413.         RGBPtr = &LP->table[theEntry].rgb;
  414.         *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  415.         _VToGo=_V - LP->_VFixed + LP->_VHalfStep;
  416.         for(i=LP->fixed;i<LP->dacs;i++) {
  417.             _g=LP->_gain[i];
  418.             Vi=_VToGo/_g;                        /* truncate to integer */
  419.             if(Vi>LP->VMax)Vi=LP->VMax;
  420.             if(Vi<LP->VMin)Vi=LP->VMin;
  421.             _VToGo -= Vi*_g;
  422.             ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  423.         }
  424.         /***************************************************************/
  425.     }
  426.     return;
  427. }
  428.  
  429. void _SetLuminance(LuminanceRecord *LPtr,int theEntry,FIXED _luminance)
  430. /*
  431. Set one entry in the ColorSpec table to a specified luminance. This is the
  432. private subroutine that actually does all the work for SetLuminance(). This
  433. routine is designed to run as fast as possible.
  434. */
  435. {
  436.     register LuminanceRecord *LP=LPtr;
  437.     register int i,Vi;
  438.     register FIXED _VToGo,_g;
  439.     RGBColor *RGBPtr;
  440.     register short leftShift=LP->leftShift;
  441.     
  442.     RGBPtr = &LP->table[theEntry].rgb;
  443.     *RGBPtr=LP->rgb;                        /* load for fixed DACs */
  444.     _VToGo=_LToV(LP,_luminance + LP->_LOffset) - LP->_VFixed + LP->_VHalfStep;
  445.     for(i=LP->fixed;i<LP->dacs;i++) {
  446.         _g=LP->_gain[i];
  447.         Vi=_VToGo/_g;                        /* truncate to integer */
  448.         if(Vi>LP->VMax)Vi=LP->VMax;
  449.         if(Vi<LP->VMin)Vi=LP->VMin;
  450.         _VToGo -= Vi*_g;
  451.         ((short *)RGBPtr)[LP->dac[i]]=Vi<<leftShift;
  452.     }
  453. }
  454.  
  455. void IncrementLuminance(GDHandle device,register LuminanceRecord *LP,int theEntry)
  456. /*
  457. Make smallest possible increase of the luminance of one entry in the ColorSpec
  458. table. If necessary, bump up the upper luminance bound to include the incremented luminance.
  459. */
  460. {
  461.     double V,L,highL;
  462.     
  463.     V=GetV(device,LP,theEntry);
  464.     V+=LP->VHalfStep;
  465.     V+=LP->VHalfStep;
  466.     L=VToL(LP,V);
  467.     if(L>LP->highLuminance)highL=L;
  468.     else highL=LP->highLuminance;
  469.     SetLuminance(device,LP,theEntry,L,LP->lowLuminance,highL);
  470. }
  471.  
  472. FIXED _Tolerance(LuminanceRecord *LPtr,FIXED _luminance)
  473. {
  474.     register LuminanceRecord *LP=LPtr;
  475.     register FIXED _tolerance,_t;
  476.     
  477.     /* quick and dirty estimate of tolerance */
  478.     if(LP->L.exists==luminanceSet){
  479.         _tolerance=LP->L._Lu[0] - _luminance;
  480.         _t=_luminance - LP->L._Lu[LUMINANCES_IN_TABLE-1];
  481.         if(_tolerance<_t)_tolerance=_t;
  482.         if(_tolerance<0)_tolerance=0;
  483.         _tolerance+=LP->_tolerance;
  484.     }
  485.     else _tolerance=LP->_tolerance;
  486.     return _tolerance;
  487. }
  488.  
  489. double GetLuminance(GDHandle device,LuminanceRecord *LP,int theEntry)
  490. /*
  491. If device is not NULL then examines one entry in the actual clut, otherwise
  492. examines the ColorSpec table contained in *LP, and in either case returns the
  493. luminance that will be produced. Supplying an illegal entry value results in a
  494. returned value of -INF.
  495. */
  496. {
  497.     double v;
  498. static LuminanceRecord *LPSave;
  499. LPSave=LP;
  500.     v=GetV(device,LP,theEntry);
  501. assert(LP==LPSave);
  502.     return VToL(LP,v);
  503. }
  504.  
  505. double GetV(GDHandle device,LuminanceRecord *LP,int theEntry)
  506. /*
  507. If device is not NULL (and has a driver) then examines one entry in the actual clut, otherwise
  508. examines the ColorSpec table contained in *LP, and in either case returns the V
  509. that will be produced. Supplying an illegal entry value results in a returned
  510. value of -INF.
  511. */
  512. {
  513.     RGBColor *myRGBPtr=NULL;
  514.     ColorSpec table[256];
  515.     double V;
  516.     int error=0;
  517.     
  518.     if(device != NULL && (*device)->gdRefNum!=0) {
  519.         if(theEntry<0 || theEntry>=GDCLUTSIZE(device)){
  520.             PrintfExit("GetLuminance: illegal entry %d\n",theEntry);
  521.         }
  522.         assert(GDCLUTSIZE(device)<=256);
  523.         error=GDGetEntries(device,0,GDCLUTSIZE(device)-1,table);
  524.         if(error){
  525.             PrintfExit("GetLuminance: GDGetEntries error %d\007",error);
  526.         }
  527.         myRGBPtr=&table[theEntry].rgb;
  528.     }else{
  529.         if(LP!=NULL && theEntry>=0 && theEntry<COLORS) 
  530.             myRGBPtr=&LP->table[theEntry].rgb;
  531.     }
  532.     if(myRGBPtr == NULL) return -INF;
  533.     V = LP->r*(myRGBPtr->red>>8);
  534.     V += LP->g*(myRGBPtr->green>>8);
  535.     V += LP->b*(myRGBPtr->blue>>8);
  536.     return V;
  537. }
  538.  
  539. void LoadLuminances(GDHandle device, LuminanceRecord *LP,
  540.     int firstEntry, int lastEntry)
  541. /*
  542. This just calls GDSetEntries() or GDDirectSetEntries() to load your ColorSpec
  543. table into the clut of your screen device. It is here simply to provide a
  544. cosmetic match to the call to SetLuminances(), for which loading the clut is
  545. optional. Note: if you prefer, instead of LP you may send just the address of a
  546. ColorSpec table, cast to (LuminanceRecord *), since a LuminanceRecord begins
  547. with a ColorSpec table.
  548. */
  549. {
  550.     int error;
  551.     short isGray,i;
  552.     RGBColor *RGBPtr;
  553.     ColorSpec table[256],*tablePtr;
  554.  
  555.     if(device==NULL || (*device)->gdRefNum==0)return;    // no device, or no device driver
  556.     isGray=!TestDeviceAttribute(device,gdDevType);
  557.     if(isGray){
  558.         for(i=firstEntry;i<=lastEntry;i++){
  559.             RGBPtr=&table[i].rgb;
  560.             RGBPtr->red=RGBPtr->green=RGBPtr->blue=LP->table[i].rgb.green;
  561.         }
  562.         tablePtr=table;
  563.     }else tablePtr=LP->table;
  564.     switch((*device)->gdType){
  565.     case fixedType:
  566.         printf("LoadLuminances: this device has a fixed CLUT.\n");
  567.         return;
  568.     case clutType:
  569.         error=GDSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr+firstEntry);
  570.         break;
  571.     case directType:
  572.         error=GDDirectSetEntries(device,firstEntry,lastEntry-firstEntry,tablePtr+firstEntry);
  573.         break;
  574.     }
  575.     if(error) printf("\007LoadLuminances: error %d\n",error);
  576. }
  577.  
  578. double LToL(LuminanceRecord *LP,double L)
  579. {
  580.     if(L > LP->LMax) return LP->LMax;
  581.     if(L < LP->LMin) return LP->LMin;
  582.     return L;
  583. }
  584.  
  585. double VToL(LuminanceRecord *LP,double V)
  586. /*
  587. Return the luminance that would result from a nominal voltage.
  588. */
  589. {
  590. assert(LP->dacSize==8);
  591.     return FixToDouble(_VToL(LP,DoubleToFix(V)));
  592. }
  593.  
  594. FIXED _VToL(LuminanceRecord *LP,FIXED _V)
  595. /*
  596. Return the luminance that would result from a nominal voltage.
  597. */
  598. {
  599.     register int i;
  600.     register LuminanceTable *LT;
  601.     register FIXED _di;
  602.     double LF,VF;
  603.     
  604.     LT=&LP->L;
  605.     if(LT->exists==luminanceSet){
  606.         _di=FixDiv(_V-LT->_VMin,LT->_dV);
  607.         i=FixToLong(_di);
  608.         _di -= LongToFix(i);
  609.         if(i<0)return LT->_Lu[0];
  610.         if(i>=LUMINANCES_IN_TABLE-1)return LT->_Lu[LUMINANCES_IN_TABLE-1];
  611.         return LT->_Lu[i]+FixMul(_di,LT->_Lu[i+1]-LT->_Lu[i]);
  612.     }
  613.     VF=FixToDouble(_V);
  614.     if(LP->powerError < 2.0*LP->polynomialError) LF=VToLPower(LP,VF);
  615.     else LF=VToLPolynomial(LP,VF);
  616.     return DoubleToFix(LF);
  617. }
  618.  
  619. double LToV(LuminanceRecord *LP,double L)
  620. {
  621.     return FixToDouble(_LToV(LP,DoubleToFix(L)));
  622. }
  623.  
  624. FIXED _LToV(LuminanceRecord *LP,FIXED _Lu)
  625. /*
  626. New, faster version uses a tabulated version of the gamma function VToL(). This
  627. replaces the old function LToVFormulaic(), formerly called LToV(). A subtlety is
  628. that I tabulate the gamma function _Lu[]=VToL() rather than its inverse, because
  629. this yields a much lower upper bound on the error of the linearly interpolated
  630. L. This is true because the gamma function is roughly parabolic. Equally spaced
  631. samples of a parabolic function yield equal maximum error of linear
  632. interpolation in all intervals. The cost of tabulating _Lu[] rather than its
  633. inverse is that we must search the table in order to find the largest _Lu[] less
  634. than or equal to L. The search time is minimized by using a search procedure
  635. that begins the search at the place found by the last search, since successive
  636. calls to LToV() are likely to request similar luminances.
  637. */
  638. {
  639.     register LuminanceTable *LT;
  640.     register int lo,hi,i;
  641.     register FIXED _dL,_dLTemp;
  642.     FIXED _V,_dLLo;
  643.     
  644.     LT=&LP->L;
  645.     if(LT->exists != luminanceSet){
  646.         /* get domain of the monotonic section of the gamma function */
  647.         LT->_VMin=DoubleToFix(LToVFormulaic(LP,LP->LMin));
  648.         LT->_VMax=DoubleToFix(LToVFormulaic(LP,LP->LMax));
  649.         LT->_dV=(LT->_VMax-LT->_VMin)/(LUMINANCES_IN_TABLE-1);
  650.         _V=LT->_VMin;
  651.         for(i=0;i<LUMINANCES_IN_TABLE;i++,_V+=LT->_dV) LT->_Lu[i]=_VToL(LP,_V);
  652.         /* just in case, impose monotonicity, from center on out */
  653.         for(i=LUMINANCES_IN_TABLE/2;i>0;i--)
  654.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i-1]=LT->_Lu[i];
  655.         for(i=1+LUMINANCES_IN_TABLE/2;i<LUMINANCES_IN_TABLE;i++)
  656.             if(LT->_Lu[i-1] > LT->_Lu[i]) LT->_Lu[i]=LT->_Lu[i-1];
  657.         #if FAST_LUMINANCE    /* compute the slope dL/dV */
  658.             /* find shift that maximizes ΔL precision without overflow */
  659.             _dL=0;
  660.             for(i=0;i<LUMINANCES_IN_TABLE-1;i++){
  661.                 _dLTemp=LT->_Lu[i+1]-LT->_Lu[i];
  662.                 if(_dL<_dLTemp)_dL=_dLTemp;
  663.             }
  664.             LT->LShift=Log2L(LONG_MAX/_dL);
  665.         #endif                
  666.         LT->exists=luminanceSet;
  667.         LT->latestIndex=-2;    /* invalid, so hunt will start from scratch */
  668.     }
  669.     /* hunt for L in LT->_Lu[] */
  670.     /* first check at latestIndex, latestIndex±1, ... */
  671.     lo=-1;
  672.     hi=LUMINANCES_IN_TABLE;
  673.     i=LT->latestIndex;
  674.     if(i>lo && i<hi){
  675.         if(_Lu>LT->_Lu[i])
  676.             {lo=i;i++;}
  677.         else
  678.             {hi=i;i--;}
  679.     }
  680.     /* simple bisection search, see Numerical Recipes in C, pages 98-99. */
  681.     while(hi-lo>1){
  682.         i=(hi+lo)>>1;
  683.         if(_Lu>LT->_Lu[i])lo=i;
  684.         else hi=i;
  685.     }
  686.     LT->latestIndex=lo;
  687.     if(lo<0){
  688.         LT->latestIndex=0;            /* nearest legal index */
  689.         return LT->_VMin;
  690.     }
  691.     if(lo>=LUMINANCES_IN_TABLE-1)return LT->_VMax;
  692.     _V=LT->_VMin+lo*LT->_dV;
  693.     _dL=_Lu-LT->_Lu[lo];
  694.     _dLLo=LT->_Lu[lo+1]-LT->_Lu[lo];
  695.     if(FAST_LUMINANCE){
  696.         if(_dL<LONG_MAX/LT->_dV)
  697.             _V+=LT->_dV*_dL/_dLLo;
  698.         else if((_dLLo<<LT->LShift) >= LT->_dV)
  699.             _V+=(_dL<<LT->LShift)/((_dLLo<<LT->LShift)/LT->_dV);
  700.     }else _V+=LT->_dV*(_dL/_dLLo);
  701.     return _V;
  702. }
  703.  
  704. double LToVFormulaic(LuminanceRecord *LP,double L)
  705. /*
  706. Find a nominal voltage that would give luminance L. Failing that, it returns the
  707. closest possible value. The answer is computed by inverting the formula for the
  708. gamma function.
  709. */
  710. {
  711.     if(LP->powerError < 2.0*LP->polynomialError) return LToVPower(LP,L);
  712.     else return LToVPolynomial(LP,L);
  713. }
  714.  
  715. double VToLPolynomial(LuminanceRecord *LP,double V)
  716. /*
  717. Return the luminance that would result from a nominal voltage. Uses a polynomial
  718. equation L = p[0] + p[1]*V^1 +... of any degree. dgp 8/11/89
  719. */
  720. {
  721.     register double L,VV;
  722.     register int i,m;
  723.     
  724.     m=LP->coefficients;
  725.     L=0.0;
  726.     VV=1.0;
  727.     for(i=0;i<m;i++){
  728.         L+=LP->p[i]*VV;
  729.         VV*=V;
  730.     }
  731.     return L;
  732. }
  733.  
  734. double VToLPower(LuminanceRecord *LP,double V)
  735. /*
  736. Return the luminance that would result from a nominal voltage.
  737. Uses a rectified power law:
  738.  
  739.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  740.     
  741.     Rectify(x)=x if x≥0
  742.     Rectify(x)=0 if x<0
  743.  
  744. dgp 10/28/89
  745. */
  746. {
  747.     double L;
  748.     
  749.     L=LP->power[1]+LP->power[2]*V;
  750.     if(L>0.0) L=LP->power[0]+pow(L,LP->power[3]);
  751.     else L=LP->power[0];
  752.     return L;
  753. }
  754.  
  755. double LToVPolynomial(LuminanceRecord *LP,double L)
  756. /*
  757. Find a nominal voltage that would give luminance L. Failing that, it returns the
  758. closest possible value. Solves a polynomial equation L = p[0] + p[1]*V^1 +... of
  759. any degree. First we find a quick approximate answer by solving the power law
  760. fit, LToVPower(). Then we use Newton's method to home in quickly on the solution
  761. to the polynomial fit. Newton's method is taken from Numerical Recipes in C. For
  762. a polynomial of degree 8 (which is what I recommend) this routine now does about
  763. 3? iterations and takes about 1.1? ms. The error in V is less than TOLERANCE,
  764. i.e. relative to full scale of 256 we get an accuracy of 14.6 bits. You can
  765. increase the TOLERANCE to reduce the number of iterations to make it slightly
  766. faster. A TOLERANCE of 1.0 (for 8-bit accuracy) takes 0.8? ms, not much of a
  767. savings in time, and a large loss in accuracy. Reducing TOLERANCE to 1e-8 yields
  768. an accuracy of nearly 35 bits and takes 1.5? ms, only slightly longer.
  769.  
  770. This routine will be called essentially once per clut entry, so computing a
  771. whole new clut with entries 0..255 will take 256*1.1? ms = 282? ms. Some
  772. experiments use only part of the clut, and will therefore take less time. In
  773. some of the experiments that do use the whole clut it may be desirable to
  774. compute the clut table in advance.
  775. dgp 8/11/89
  776. */
  777. {
  778.     #define TOLERANCE 0.01                    /* desired accuracy of V, see above */
  779.     #define JMAX 30
  780.     double a[MAX_COEFFICIENTS];
  781.     register int i;
  782.     int m,j;
  783.     register double f,df,u,VV,V,dV;
  784.     
  785.     if(L>LP->LMax) return LP->VMax;
  786.     if(L<LP->LMin) return LP->VMin;
  787.  
  788.     m=LP->coefficients;
  789.     if(m>MAX_COEFFICIENTS || m<1)
  790.         PrintfExit("LToVPolynomial: %d coefficients is too many or too few\007\n",m);
  791.     for(i=0;i<m;i++) a[i]=LP->p[i];
  792.     a[0]-=L;
  793.     if(LP->powerError < 0.2*LP->LMax) V=LToVPower(LP,L);    /* very good initial guess */
  794.     else {
  795.         if(LP->quadraticError < 0.2*LP->LMax) V=LToVQuadratic(LP,L);/* fair initial guess */
  796.         else PrintfExit("LToVPolynomial: neither power nor quadratic fit is available\007\n");
  797.     }
  798.     for(j=0;j<JMAX;j++){
  799.         /* evaluate function, i.e. the polynomial, and its derivative */
  800.         f=a[0];
  801.         df=0.0;
  802.         VV=1.0;
  803.         for(i=1;i<m;i++){
  804.             u=a[i]*VV;
  805.             df+=i*u;
  806.             f+=V*u;
  807.             VV*=V;
  808.         }
  809.         dV=f/df;    /* apply Newton's method */
  810.         if(V-dV<LP->VMin) dV=(V-LP->VMin)/2.0;    /* bisect */
  811.         if(V-dV>LP->VMax) dV=(V-LP->VMax)/2.0;    /* bisect */
  812.         V -= dV;
  813.         if(fabs(dV) < TOLERANCE) return V;
  814.     }
  815.     printf("LToVPolynomial(%7.2f): Warning, too many iterations. VToL(%5.1f)=%7.2f\n",V,L,VToL(LP,V));
  816.     return V;
  817.     #undef TOLERANCE
  818.     #undef JMAX
  819. }
  820.  
  821. double LToVPower(LuminanceRecord *LP,double L)
  822. /*
  823. Find a nominal voltage that would give luminance L. Failing that, it returns the
  824. closest possible value. Solves the rectified power law:
  825.     L(V)=power[0]+Rectify(power[1]+power[2]*V)^power[3]
  826.     Rectify(x)=x if x≥0
  827.     Rectify(x)=0 if x<0
  828. LToVPower takes about 300 microseconds on a Mac II with 8881 arithmetic. This
  829. routine will be called essentially once per clut entry, so computing a whole new
  830. clut with entries 0..255 will take 256*0.3 ms = 77 ms. Some experiments use only
  831. part of the clut, and will therefore take less time. In experiments that do use
  832. the whole clut it may be desirable to compute the clut table in advance.
  833. dgp 10/28/89
  834. */
  835. {
  836.     double V;
  837.     
  838.     if(L<LP->power[0]) L=LP->power[0];
  839.     V=(pow(L-LP->power[0],1.0/LP->power[3])-LP->power[1])/LP->power[2];
  840.     if(V>LP->VMax) V=LP->VMax;
  841.     return V;
  842. }
  843.  
  844. double LToVQuadratic(LuminanceRecord *LP,double L)
  845. /*
  846. This is a quick, approximate version of LToV() that is used by LToVPolynomial
  847. for an initial guess if no power law fit is available. Find nominal voltage that
  848. would give luminance L. Since a quadratic fit to the luminance calibration is
  849. only a fair approximation, this routine serves only to find a quick approximate
  850. answer to the root of a higher order polynomial fit. Solves the quadratic
  851. equation by method recommended in Numerical Recipes in C. My choice of root
  852. assumes that L is an increasing function of V over the operating range of the
  853. display. If concave (i.e. q[2]<0) use the larger root. If convex (i.e. q[2]>0)
  854. use the smaller root. In plain English, if the parabola is right side up then we
  855. take the right hand branch. If the parabola is upside down then we take the left
  856. hand branch. dgp 8/12/89
  857. */
  858. {
  859.     register double a,b,c,q;
  860.     double V;
  861.  
  862.     if(L>LP->LMax) return LP->VMax;
  863.     if(L<LP->LMin) return LP->VMin;
  864.     c=LP->q[0]-L;
  865.     b=LP->q[1];
  866.     a=LP->q[2];
  867.     if(a == 0.0) return -c/b;
  868.     if(b<0.0)
  869.         q=-0.5*(b-sqrt(b*b-4.0*a*c));
  870.     else
  871.         q=-0.5*(b+sqrt(b*b-4.0*a*c));
  872.     if(a*b<0.0)
  873.         V=q/a;
  874.     else
  875.         V=c/q;
  876.     return V;
  877. }
  878.  
  879. double SetLuminanceRange(LuminanceRecord *LPtr,double lowLuminance,double highLuminance)
  880. /*
  881. The user will never call this routine directly, as it is called automatically by
  882. SetLuminance and SetLuminances. This routine does the hard work of figuring out
  883. which DACs to fix to optimally represent a given luminance range. The goal is to
  884. fix the coarser DACs and vary the finer DACs to represent the image. The CHANGES
  885. in luminance (i.e. the contrast) will have the accuracy of the coarsest of the
  886. variable DACs. To avoid the overhead of figuring all this out every single time
  887. you call SetLuminance(), the results are stored temporarily in your
  888. LuminanceRecord, and are only recomputed if you request a different luminance
  889. range.
  890.  
  891. This function checks and returns immediately if you call it with the same
  892. luminance range that it already has stored in your LuminanceRecord.
  893. */
  894. {
  895.     register LuminanceRecord *LP=LPtr;
  896.     register int i;
  897.     double dL,dL2,LOffset;
  898.     double VSum,VFixed;
  899.     int dacs,fixed,dac[DACS];
  900.     double gain[DACS],V[DACS];
  901.     double VGoal;
  902.     double lowV,highV;
  903.     double tolerance;
  904.     double VHalfStep;
  905.     double dV;
  906.     unsigned short *RGBArray;
  907.     
  908.     /* -1. If necessary, swap low & high */
  909.     if(lowLuminance > highLuminance){
  910.         dL=lowLuminance;
  911.         lowLuminance=highLuminance;
  912.         highLuminance=dL;
  913.     }
  914.  
  915.     /* 0. Return right away if our work has already been done */
  916.     if(LP->rangeSet == luminanceSet &&
  917.         LP->lowLuminance==lowLuminance &&
  918.         LP->highLuminance==highLuminance)
  919.             return LP->tolerance;
  920.  
  921.     /* 1. sort the DACs by gain, discard any with zero or negative gain.
  922.     Gain is defined as the change in V when a single DAC is increased from 0 to 1.
  923.     Sort the gains so that gain[0]≥gain[1]≥gain[2]...
  924.     Note this order is opposite to that used in Pelli & Zhang (1991). The answers
  925.     are the same, but the notation is different.
  926.     */
  927.     if(LP->rangeSet != luminanceSet){
  928.         /* This only needs to be done once, even if the range changes. */
  929.         gain[0]=LP->r;
  930.         gain[1]=LP->g;
  931.         gain[2]=LP->b;
  932.         #if DACS>3
  933.             #error "Current implementation doesn't support more than 3 DACs."
  934.         #endif
  935.         for(i=0;i<DACS;i++)dac[i]=i;
  936.         Sort(DACS,gain,dac);            /* sort so that gain[0]≥gain[1]≥gain[2] */
  937.         dacs=0;
  938.         for(i=0;i<DACS;i++) if(gain[i]>0.0)dacs++;
  939.         VHalfStep=gain[dacs-1]/2.0;        /* half the smallest step */
  940.         for(i=0;i<COLORS;i++)LP->table[i].value=i;
  941.         LP->dacSize=Log2L(LP->VMax*2-1);    // in bits, rounded up
  942.         if(LP->dacSize!=8)printf(
  943.             "Caution: Your video card seems to have %ld-bit video dacs, and Luminance.c\n"
  944.             "has only been tested with 8-bit dacs.\n",LP->dacSize);
  945.         /* Apple's convention is to provide a 16-bit value to the dac. Luminance.c
  946.         saves time and space by computing values with only as much precision as
  947.         needed for the the video card's dac (i.e. 0 to VMax, or dacSize bits). This
  948.         value must be shifted left by leftShift bits to yield a standard 16-bit value.
  949.         Note that Apple normally recommends scaling the value to fill the entire
  950.         16-bit range, e.g. multiplying an 8-bit value by 0x101, but that doesn't
  951.         seem appropriate here, where we know the dac size, and are more concerned with
  952.         step size that with range.
  953.         */
  954.         LP->leftShift=16-LP->dacSize;
  955.     }
  956.     else {
  957.         for(i=0;i<DACS;i++){
  958.             gain[i]=LP->gain[i];
  959.             dac[i]=LP->dac[i];
  960.         }
  961.         dacs=LP->dacs;
  962.         VHalfStep=LP->VHalfStep;
  963.     }
  964.     for(i=0;i<DACS;i++)V[i]=0.0;    /* initialize fixed dac voltages */
  965.     
  966.     /* 2. Transform lowLuminance and highLuminance to lowV and highV. */
  967.     lowV=LToV(LP,lowLuminance);        /* takes 0.15 ms */
  968.     highV=LToV(LP,highLuminance);    /* takes 0.15 ms */
  969.  
  970.     /* 3. Designate the finest dacs as variable until we have enough to cover
  971.     the range, then designate the rest as fixed. (DACs 0..fixed-1 will be fixed.)
  972.     */
  973.     VGoal=fabs(highV-lowV);
  974.     VSum=0.0;
  975.     fixed=0;
  976.     for(i=dacs-1;i>=0;i--) {
  977.         if(VSum >= VGoal)fixed++;
  978.         VSum += gain[i]*(LP->VMax - LP->VMin);
  979.     }
  980.     /* sum gains of all variable dacs. This is called gVary by Pelli & Zhang (1991) */
  981.     tolerance=0.0;
  982.     for(i=dacs-1;i>=fixed;i--)tolerance+=gain[i];
  983.     /* scale by highest luminance gain in the requested range */
  984.     dL=fabs(VToL(LP,lowV+1.0)-VToL(LP,lowV));
  985.     dL2=fabs(VToL(LP,highV)-VToL(LP,highV-1.0));
  986.     tolerance *= MAX(dL,dL2);
  987.  
  988.     /* 4. Temporarily set the variable DACs to their midpoints. Now set the fixed DACs
  989.     to most accurately represent (lowV+highV)/2.
  990.     */
  991.     VGoal=(lowV+highV)/2.0;
  992.     VSum=0.0;
  993.     VFixed=(LP->VMax + LP->VMin)/2.0;
  994.     for(i=fixed;i<dacs;i++) VSum+=gain[i]*VFixed;    /* The mid-voltage of the variable DACs */
  995.     for(i=0;i<fixed;i++) {
  996.         V[i]=floor((VHalfStep+VGoal-VSum)/gain[i]);    /* the unattenuated voltage of DAC i */
  997.         V[i]=MAX(LP->VMin,MIN(LP->VMax,V[i]));
  998.         VSum += V[i]*gain[i];
  999.     }
  1000.     VSum=0.0;
  1001.     for(i=0;i<fixed;i++) VSum+=gain[i]*V[i];
  1002.     VFixed=VSum;
  1003.     
  1004.     /* 5. The limited precision of the fixed DACs may result in a 
  1005.     small offset between the requested and now available range. The offset will be
  1006.     at most half a step of the finest of the fixed DACs. It seems reasonable
  1007.     to offset the requested luminance range by up to that small amount to better fit
  1008.     the now available range. */
  1009.     LOffset=0.0;
  1010.     if(fixed>0){
  1011.         VSum=VFixed;
  1012.         for(i=fixed;i<dacs;i++) VSum += LP->VMin*gain[i];
  1013.         dV= VSum - MIN(lowV,highV);
  1014.         dV= MIN(dV,gain[fixed-1]/2.0);
  1015.         if(dV>0.0) LOffset+=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1016.         VSum=VFixed;
  1017.         for(i=fixed;i<dacs;i++) VSum += LP->VMax*gain[i];
  1018.         dV= MAX(lowV,highV) - VSum;
  1019.         dV= MIN(dV,gain[fixed-1]/2.0);
  1020.         if(dV>0.0) LOffset-=VToL(LP,VSum+dV)-VToL(LP,VSum);
  1021.     }
  1022.     /* Now save this information in the LuminanceRecord for future use */
  1023.     for(i=0;i<DACS;i++){
  1024.         LP->dac[i]=dac[i];
  1025.         LP->gain[i]=gain[i];
  1026.     }
  1027.     LP->rangeSet=luminanceSet;
  1028.     LP->lowLuminance=lowLuminance;
  1029.     LP->highLuminance=highLuminance;
  1030.     LP->dacs=dacs;
  1031.     LP->VHalfStep=VHalfStep;
  1032.     LP->fixed=fixed;
  1033.     LP->VFixed=VFixed;
  1034.     LP->tolerance=tolerance;
  1035.     LP->LOffset=LOffset;
  1036.  
  1037.     /* Copy into FIXED variables */
  1038.     for(i=0;i<DACS;i++)LP->_gain[i]=DoubleToFix(LP->gain[i]);
  1039.     LP->_VHalfStep=DoubleToFix(VHalfStep);
  1040.     LP->_VFixed=DoubleToFix(VFixed);
  1041.     LP->_tolerance=DoubleToFix(tolerance);
  1042.     LP->_LOffset=DoubleToFix(LOffset);
  1043.  
  1044.     /* cache the fixed DAC values, and zero the rest for good measure */
  1045.     RGBArray = (unsigned short *) &LP->rgb;    /* treat as an array */
  1046.     for(i=0;i<DACS;i++){
  1047.         RGBArray[LP->dac[i]] = (short)V[i]<<LP->leftShift;
  1048.     }
  1049.     return tolerance;
  1050. }
  1051.  
  1052. static void Sort(int n,double arr[],int krr[])
  1053. /*
  1054. Slightly modified sort routine piksr2() from Numerical Recipes in C. I changed
  1055. "float" to "double", and I changed second array to int. I reversed the order, so
  1056. largest element would be first. I changed it to use conventional c arrays,
  1057. starting at index 0.
  1058. */
  1059. {
  1060.     register int i,j;
  1061.     double a;
  1062.     int k;
  1063.  
  1064.     for(j=1;j<n;j++) {
  1065.         a=arr[j];
  1066.         k=krr[j];
  1067.         i=j-1;
  1068.         while (i >= 0 && arr[i] < a) {
  1069.             arr[i+1]=arr[i];
  1070.             krr[i+1]=krr[i];
  1071.             i--;
  1072.         }
  1073.         arr[i+1]=a;
  1074.         krr[i+1]=k;
  1075.     }
  1076. }
  1077.